COMMENT ⊗   VALID 00002 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	      SUBROUTINE FMFP(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H)
C00016 ENDMK
C⊗;
      SUBROUTINE FMFP(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H)
C
C     ..................................................................
C
C	 SUBROUTINE FMFP
C
C	 PURPOSE
C	    TO FIND A LOCAL MINIMUM OF A FUNCTION OF SEVERAL VARIABLES
C	    BY THE METHOD OF FLETCHER AND POWELL
C
C	 USAGE
C	    CALL FMFP(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H)
C
C	 DESCRIPTION OF	PARAMETERS
C	    FUNCT  - USER-WRITTEN SUBROUTINE CONCERNING	THE FUNCTION TO
C		     BE	MINIMIZED. IT MUST BE OF THE FORM
C		     SUBROUTINE	FUNCT(N,ARG,VAL,GRAD)
C		     AND MUST SERVE THE	FOLLOWING PURPOSE
C		     FOR EACH N-DIMENSIONAL ARGUMENT VECTOR  ARG,
C		     FUNCTION VALUE AND	GRADIENT VECTOR	MUST BE	COMPUTED
C		     AND, ON RETURN, STORED IN VAL AND GRAD RESPECTIVELY
C	    N	   - NUMBER OF VARIABLES
C	    X	   - VECTOR OF DIMENSION N CONTAINING THE INITIAL
C		     ARGUMENT WHERE THE	ITERATION STARTS. ON RETURN,
C		     X HOLDS THE ARGUMENT CORRESPONDING	TO THE
C		     COMPUTED MINIMUM FUNCTION VALUE
C	    F	   - SINGLE VARIABLE CONTAINING	THE MINIMUM FUNCTION
C		     VALUE ON RETURN, I.E. F=F(X).
C	    G	   - VECTOR OF DIMENSION N CONTAINING THE GRADIENT
C		     VECTOR CORRESPONDING TO THE MINIMUM ON RETURN,
C		     I.E. G=G(X).
C	    EST	   - IS	AN ESTIMATE OF THE MINIMUM FUNCTION VALUE.
C	    EPS	   - TESTVALUE REPRESENTING THE	EXPECTED ABSOLUTE ERROR.
C		     A REASONABLE CHOICE IS 10**(-6), I.E.
C		     SOMEWHAT GREATER THAN 10**(-D), WHERE D IS	THE
C		     NUMBER OF SIGNIFICANT DIGITS IN FLOATING POINT
C		     REPRESENTATION.
C	    LIMIT  - MAXIMUM NUMBER OF ITERATIONS.
C	    IER	   - ERROR PARAMETER
C		     IER = 0 MEANS CONVERGENCE WAS OBTAINED
C		     IER = 1 MEANS NO CONVERGENCE IN LIMIT ITERATIONS
C		     IER =-1 MEANS ERRORS IN GRADIENT CALCULATION
C		     IER = 2 MEANS LINEAR SEARCH TECHNIQUE INDICATES
C		     IT	IS LIKELY THAT THERE EXISTS NO MINIMUM.
C	    H	   - WORKING STORAGE OF	DIMENSION N*(N+7)/2.
C
C	 REMARKS
C	     I)	THE SUBROUTINE NAME REPLACING THE DUMMY	ARGUMENT  FUNCT
C		MUST BE	DECLARED AS EXTERNAL IN	THE CALLING PROGRAM.
C	    II)	IER IS SET TO 2	IF , STEPPING IN ONE OF	THE COMPUTED
C		DIRECTIONS, THE	FUNCTION WILL NEVER INCREASE WITHIN
C		A TOLERABLE RANGE OF ARGUMENT.
C		IER = 2	MAY OCCUR ALSO IF THE INTERVAL WHERE F
C		INCREASES IS SMALL AND THE INITIAL ARGUMENT WAS
C		RELATIVELY FAR AWAY FROM THE MINIMUM SUCH THAT THE
C		MINIMUM	WAS OVERLEAPED.	THIS IS	DUE TO THE SEARCH
C		TECHNIQUE WHICH	DOUBLES	THE STEPSIZE UNTIL A POINT
C		IS FOUND WHERE THE FUNCTION INCREASES.
C
C	 SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C	    FUNCT
C
C	 METHOD
C	    THE	METHOD IS DESCRIBED IN THE FOLLOWING ARTICLE
C	    R. FLETCHER	AND M.J.D. POWELL, A RAPID DESCENT METHOD FOR
C	    MINIMIZATION,
C	    COMPUTER JOURNAL VOL.6, ISS. 2, 1963, PP.163-168.
C
C     ..................................................................
C
C
C	 DIMENSIONED DUMMY VARIABLES
      DIMENSION	H(1),X(1),G(1)
C
C	 COMPUTE FUNCTION VALUE	AND GRADIENT VECTOR FOR	INITIAL	ARGUMENT
      CALL FUNCT(N,X,F,G)
C
C	 RESET ITERATION COUNTER AND GENERATE IDENTITY MATRIX
      IER=0
      KOUNT=0
      N2=N+N
      N3=N2+N
      N31=N3+1
    1 K=N31
      DO 4 J=1,N
      H(K)=1.
      NJ=N-J
      IF(NJ)5,5,2
    2 DO 3 L=1,NJ
      KL=K+L
    3 H(KL)=0.
    4 K=KL+1
C
C	 START ITERATION LOOP
    5 KOUNT=KOUNT +1
C
C	 SAVE FUNCTION VALUE, ARGUMENT VECTOR AND GRADIENT VECTOR
      OLDF=F
      DO 9 J=1,N
      K=N+J
      H(K)=G(J)
      K=K+N
      H(K)=X(J)
C
C	 DETERMINE DIRECTION VECTOR H
      K=J+N3
      T=0.
      DO 8 L=1,N
      T=T-G(L)*H(K)
      IF(L-J)6,7,7
    6 K=K+N-L
      GO TO 8
    7 K=K+1
    8 CONTINUE
    9 H(J)=T
C
C	 CHECK WHETHER FUNCTION	WILL DECREASE STEPPING ALONG H.
      DY=0.
      HNRM=0.
      GNRM=0.
C
C	 CALCULATE DIRECTIONAL DERIVATIVE AND TESTVALUES FOR DIRECTION
C	 VECTOR	H AND GRADIENT VECTOR G.
      DO 10 J=1,N
      HNRM=HNRM+ABS(H(J))
      GNRM=GNRM+ABS(G(J))
   10 DY=DY+H(J)*G(J)
C
C	 REPEAT	SEARCH IN DIRECTION OF STEEPEST	DESCENT	IF DIRECTIONAL
C	 DERIVATIVE APPEARS TO BE POSITIVE OR ZERO.
      IF(DY)11,51,51
C
C	 REPEAT	SEARCH IN DIRECTION OF STEEPEST	DESCENT	IF DIRECTION
C	 VECTOR	H IS SMALL COMPARED TO GRADIENT	VECTOR G.
   11 IF(HNRM/GNRM-EPS)51,51,12
C
C	 SEARCH	MINIMUM	ALONG DIRECTION	H
C
C	 SEARCH	ALONG H	FOR POSITIVE DIRECTIONAL DERIVATIVE
   12 FY=F
      ALFA=2.*(EST-F)/DY
      AMBDA=1.
C
C	 USE ESTIMATE FOR STEPSIZE ONLY	IF IT IS POSITIVE AND LESS THAN
C	 1. OTHERWISE TAKE 1. AS STEPSIZE
      IF(ALFA)15,15,13
   13 IF(ALFA-AMBDA)14,15,15
   14 AMBDA=ALFA
   15 ALFA=0.
C
C	 SAVE FUNCTION AND DERIVATIVE VALUES FOR OLD ARGUMENT
   16 FX=FY
      DX=DY
C
C	 STEP ARGUMENT ALONG H
      DO 17 I=1,N
   17 X(I)=X(I)+AMBDA*H(I)
C
C	 COMPUTE FUNCTION VALUE	AND GRADIENT FOR NEW ARGUMENT
      CALL FUNCT(N,X,F,G)
      FY=F
C
C	 COMPUTE DIRECTIONAL DERIVATIVE	DY FOR NEW ARGUMENT.  TERMINATE
C	 SEARCH, IF DY IS POSITIVE. IF DY IS ZERO THE MINIMUM IS FOUND
      DY=0.
      DO 18 I=1,N
   18 DY=DY+G(I)*H(I)
      IF(DY)19,36,22
C
C	 TERMINATE SEARCH ALSO IF THE FUNCTION VALUE INDICATES THAT
C	 A MINIMUM HAS BEEN PASSED
   19 IF(FY-FX)20,22,22
C
C	 REPEAT	SEARCH AND DOUBLE STEPSIZE FOR FURTHER SEARCHES
   20 AMBDA=AMBDA+ALFA
      ALFA=AMBDA
C	 END OF	SEARCH LOOP
C
C	 TERMINATE IF THE CHANGE IN ARGUMENT GETS VERY LARGE
      IF(HNRM*AMBDA-1.E10)16,16,21
C
C	 LINEAR	SEARCH TECHNIQUE INDICATES THAT	NO MINIMUM EXISTS
   21 IER=2
      RETURN
C
C	 INTERPOLATE CUBICALLY IN THE INTERVAL DEFINED BY THE SEARCH
C	 ABOVE AND COMPUTE THE ARGUMENT	X FOR WHICH THE	INTERPOLATION
C	 POLYNOMIAL IS MINIMIZED
   22 T=0.
   23 IF(AMBDA)24,36,24
   24 Z=3.*(FX-FY)/AMBDA+DX+DY
      ALFA=AMAX1(ABS(Z),ABS(DX),ABS(DY))
      DALFA=Z/ALFA
      DALFA=DALFA*DALFA-DX/ALFA*DY/ALFA
      IF(DALFA)51,25,25
   25 W=ALFA*SQRT(DALFA)
      ALFA=DY-DX+W+W
      IF(ALFA) 250,251,250
  250 ALFA=(DY-Z+W)/ALFA
      GO TO 252
  251 ALFA=(Z+DY-W)/(Z+DX+Z+DY)
  252 ALFA=ALFA*AMBDA
      DO 26 I=1,N
   26 X(I)=X(I)+(T-ALFA)*H(I)
C
C	 TERMINATE, IF THE VALUE OF THE	ACTUAL FUNCTION	AT X IS	LESS
C	 THAN THE FUNCTION VALUES AT THE INTERVAL ENDS.	OTHERWISE REDUCE
C	 THE INTERVAL BY CHOOSING ONE END-POINT	EQUAL TO X AND REPEAT
C	 THE INTERPOLATION.  WHICH END-POINT IS	CHOOSEN	DEPENDS	ON THE
C	 VALUE OF THE FUNCTION AND ITS GRADIENT	AT X
C
      CALL FUNCT(N,X,F,G)
      IF(F-FX)27,27,28
   27 IF(F-FY)36,36,28
   28 DALFA=0.
      DO 29 I=1,N
   29 DALFA=DALFA+G(I)*H(I)
      IF(DALFA)30,33,33
   30 IF(F-FX)32,31,33
   31 IF(DX-DALFA)32,36,32
   32 FX=F
      DX=DALFA
      T=ALFA
      AMBDA=ALFA
      GO TO 23
   33 IF(FY-F)35,34,35
   34 IF(DY-DALFA)35,36,35
   35 FY=F
      DY=DALFA
      AMBDA=AMBDA-ALFA
      GO TO 22
C
C	 TERMINATE, IF FUNCTION	HAS NOT	DECREASED DURING LAST ITERATION
   36 IF(OLDF-F+EPS)51,38,38
C
C	 COMPUTE DIFFERENCE VECTORS OF ARGUMENT	AND GRADIENT FROM
C	 TWO CONSECUTIVE ITERATIONS
   38 DO 37 J=1,N
      K=N+J
      H(K)=G(J)-H(K)
      K=N+K
   37 H(K)=X(J)-H(K)
C
C	 TEST LENGTH OF	ARGUMENT DIFFERENCE VECTOR AND DIRECTION VECTOR
C	 IF AT LEAST N ITERATIONS HAVE BEEN EXECUTED. TERMINATE, IF
C	 BOTH ARE LESS THAN  EPS
      IER=0
      IF(KOUNT-N)42,39,39
   39 T=0.
      Z=0.
      DO 40 J=1,N
      K=N+J
      W=H(K)
      K=K+N
      T=T+ABS(H(K))
   40 Z=Z+W*H(K)
      IF(HNRM-EPS)41,41,42
   41 IF(T-EPS)56,56,42
C
C	 TERMINATE, IF NUMBER OF ITERATIONS WOULD EXCEED  LIMIT
   42 IF(KOUNT-LIMIT)43,50,50
C
C	 PREPARE UPDATING OF MATRIX H
   43 ALFA=0.
      DO 47 J=1,N
      K=J+N3
      W=0.
      DO 46 L=1,N
      KL=N+L
      W=W+H(KL)*H(K)
      IF(L-J)44,45,45
   44 K=K+N-L
      GO TO 46
   45 K=K+1
   46 CONTINUE
      K=N+J
      ALFA=ALFA+W*H(K)
   47 H(J)=W
C
C	 REPEAT	SEARCH IN DIRECTION OF STEEPEST	DESCENT	IF RESULTS
C	 ARE NOT SATISFACTORY
      IF(Z*ALFA)48,1,48
C
C	 UPDATE	MATRIX H
   48 K=N31
      DO 49 L=1,N
      KL=N2+L
      DO 49 J=L,N
      NJ=N2+J
      H(K)=H(K)+H(KL)*H(NJ)/Z-H(L)*H(J)/ALFA
   49 K=K+1
      GO TO 5
C	 END OF	ITERATION LOOP
C
C	 NO CONVERGENCE	AFTER  LIMIT  ITERATIONS
   50 IER=1
      RETURN
C
C	 RESTORE OLD VALUES OF FUNCTION	AND ARGUMENTS
   51 DO 52 J=1,N
      K=N2+J
   52 X(J)=H(K)
      CALL FUNCT(N,X,F,G)
C
C	 REPEAT	SEARCH IN DIRECTION OF STEEPEST	DESCENT	IF DERIVATIVE
C	 FAILS TO BE SUFFICIENTLY SMALL
      IF(GNRM-EPS)55,55,53
C
C	 TEST FOR REPEATED FAILURE OF ITERATION
   53 IF(IER)56,54,54
   54 IER=-1
      GOTO 1
   55 IER=0
   56 RETURN
      END